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ABSTRACT 


Image  theory  is  an  ideal  method  for  calculating  the  transmission  loss  in  a  shallow 
water  (wedge  shaped  ocean)  environment.  It  can  be  used  in  cross-slope,  at  all 
frequencies  and  in  transitional  cut  off  regions  that  are  out  of  bounds  to  normal  mode 
theories. 

This  thesis  had  three  objectives;  1)  convert  the  existing  image  theory  models  called 
URTEXT  and  WEDGE  into  a  high  level  scripting  language  called  MATLAB™  by 
Math  Works,  2)  linearize  the  existing  quadruplet  expansion  program  to  increase  speed, 
and  3)  to  incorporate  the  Arctan  approximation  of  the  Rayleigh  reflection  coefficient 
into  the  quadruplet  expansion  for  the  fast  bottom  case. 

Objective  1  was  conq)leted  with  accurate  results.  Objective  2  was  completed  with 
a  factor  of  8  increase  in  speed.  Objective  3  incorporated  the  Arctan  approximation  of 
the  reflection  coefficient  for  a  fast  bottom  into  the  quadruplet  expansion,  but  due  to  the 
inaccuracy  of  the  reflection  coefficient  after  the  second  quadruplet,  the  results  were  not 
favorable.  It  was  also  discovered  that  even  widi  the  Rayleigh  reflection  coefficient,  the 
first  order  iqiptoximations  made  in  developing  the  quadruplet  expansion  equation 
(Equation  6-27)  are  not  accurate  enough  for  the  fast  bottom  case. 
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I.  INTRODUCTION 


The  current  emphasis  of  undersea  warfare  is  on  the  shallow  coastal  regions  in  littoral 
wjUers.  The  sound  propagation  in  these  regions  is  influenced  by  bottom  and  surface 
interaction  and  not  by  the  convergence  zones  and  ducting  found  in  deep  water.  The 
variables  in  shallow  water  are  bottom  type,  bottom  density,  slope,  sea  state,  surface  winds, 
sea  surface  ten^rature,  etc.,  all  of  which  can  change  quickly  or  be  unavailable  in  the  area 
of  interest. 

There  is  no  closed  form  analytic  solution  to  the  problem  of  a  penetrable,  sloping 
bottom  as  found  on  the  continental  shelf  region  or  shoreline  areas.  There  are  computer 
models  using  various  approximation  methods  ranging  from  purely  numerical  (parabolic 
equation)  to  approximately  analytical  (adiabatic  mode  theory).  Transmission  loss 
modeling  in  a  wedge  shaped  ocean  falls  into  three  major  areas:  1)  parabolic  approximation, 
2)  adiabatic  normal  mode  theory,  and  3)  image  theory. 

A.  PARABOLIC  EQUATION 

The  parabolic  equation  (PE)  is  a  range-dependent,  underwater  acoustic  propagation 
model.  In  it  the  Helmholz  equation  is  replaced  by  a  one-way  parabolic  equation  which 
generates  an  acoustic  field  as  an  initial  value  {noblem.  The  PE  ^jproach  is  limited  by 
excessive  computer  tunning  time  in  shallow  water  environments,  higher  frequencies,  and 
horizontal  rays  <  40*. 

The  first  plication  of  PE  to  underwater  acoustics  was  by  Tappert  in  1977  with  a 
restriction  on  the  maximum  angular  iqierture  of  ±20'at  the  source  [Ref.  1].  This  restriction 
was  relaxed  to  ±40*  with  a  higher  order  parabolic  equation  by  Collins  [Ref.  2]. 
Computerized  calculations  by  Jensen  and  Kuperman  [Ref.  3]  based  on  a  split-step  Fast 
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Fourier  Transform  (FFT)  solution  of  the  parabolic  equation  provided  detailed  information 
about  the  beams  projected  into  the  bottom,  showing  modal  cutoffs.  The  results  compared 
favorably  to  model  tank  experiments  by  Coppens  and  Sanders  [Ref  4].  Since  the  PE  is 
cylindrically  symmetric,  sloping  bottoms  cause  a  problem.  Lee  and  McDaniel  [Ref  5] 
sectioned  the  wedge  into  a  series  of  range  independent  regions  and  then  applied  normal  PE 
methods.  Another  approach  to  this  problem  is  by  Collins  [Ref  6]  with  the  rotated  PE. 

The  rotated  PE  steps  parallel  to  the  ocean  bottom  and  has  two  normal  derivatives 
preserving  pressure  and  particle  velocity  normal  components  thus  handling  sloping 
interfaces  properly.  Fawcett  [Ref  7]  developed  a  three-dimensional  computer  model 
using  the  wide-angle  PE  approximation  by  Thomson  [Ref  8]  with  FFT  algorithm  as  an 
azimuthal  operator.  Most  recently,  Collins  [Ref.  9]  has  developed  the  energy  conserving 
PE  for  elastic  media  to  improve  accuracy  in  range-dependent  elastic  media.  The  energy- 
conserving  elastic  PE  is  a  generalization  of  the  energy  conserving  acoustic  PE  [Ref  10].  It 
involves  approximating  a  range-dependent  waveguide  into  a  sequence  of  range- 
independent  regions  with  a  linear  iq)proximation  of  compressional  energy  flux  between  the 
vertical  interfaces. 

B.  NORMAL  MODE  THEORY 

Normal  mode  theory  provides  an  exact  solution  to  the  wave  equation.  The  distinct 
advantage  of  normal  modes  is  that  once  the  set  of  eigenfunctions  have  been  determined,  the 
range  and  depth  dependence  of  the  transmission  loss  can  be  calculated  directly  [Rwf  1 1]. 

Normal  mode  theory  is  a  range-independent  approach  but  for  a  wedge  sh^)ed  ocean,  a 
range-dependent  a{q»oach  is  needed  due  to  the  bottom  interaction  (in  the  case  of  a  lossy 
bottom)  and  bottom  reflection  angle.  Pierce  [Ref  12]  used  adiabatic  separaticm  of  depth 
and  range  in  the  wave  equation  to  get  an  approximation.  Graves,  Nagl,  Uberall,  and  Zarur 
[Ref  13]  applied  this  method  to  the  wedge  shaped  ocean  using  rigid  bottom  and  isospeed 
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water.  This  method  was  good  only  for  small  slopes.  Buckingham  [Ref.  14]  developed  a 
solution  for  the  penetrable  bottom  using  an  "effective"  wedge  with  pressure  release  bottom 
below  the  actual  penetrable  bottom.  By  using  a  coordinate  transformation,  range- 
dependent  modes  are  matched  to  range  independent  eigenfunctions  in  the  "effective"  wedge 
to  account  for  amplitude  and  phase  shifts.  Then  a  sum  of  range  dependent  normal  modes 
are  applied  to  the  new  bottom. 

A  problem  with  adiabatic  normal  mode  approximation  is  it  does  not  explain  the 
transition  to  evanescent  modes  at  cutoff.  Pierce  [Ref.  15]  combined  PE  with  adiabatic 
normal  mode  theory  to  develop  a  critical  depth  function.  Known  as  the  augmented 
adiabatic  mode  theory,  it  is  a  sum  of  modal  terms  with  each  mode  eventually  encountering 
a  critical  depth  at  a  critical  range.  The  results  compare  well  with  the  results  of  Jensen  and 
Kuperman  [Ref.  3]. 

Another  method  is  to  combine  local  modes  with  ray  acoustics.  Arnold  and  Felsen 
[Ref.  16]  develq)ed  intrinsic  wave  functions  that  are  determined  by  waveguide  geometry. 
Modes  are  mutually  decoupled  thus  propt^ate  independent  of  source  geometry  (physical 
shape)  when  superpositioned. 

C.  IMAGE  THEORY 

Image  theory  in  acoustics  is  similar  to  image  theory  in  optics.  A  source  radiating  in 
the  water  has  a  virtual  image  reflected  equidistant  above  the  water  180°  out  of  phase.  In  an 
isospeed,  wedge  sh^ted  ocean  numerous  images  are  produced  in  kaleidoscope  fashion. 

In  1966,  MaqdierstHi  and  Dainteth  [Ref.  17]  proposed  a  phase  incoherent  model  for 
upslope  i»x>pagation .  Kawamura  and  loarmou  [Ref.  18]  predicted  the  pressure  amplitude 
and  i^iase  of  the  sound  field  along  the  bottom  of  a  wedge  shrq)ed  fluid  layer  overlapping  a 
fast  fluid  bottom.  Coppens,  Humphries,  and  Sanders  [Ref.  19]  developed  a  phase  coherent 
model  to  calculate  the  pressure  amplitude  along  the  bottom  of  a  wedge  of  water  overlying  a 
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fast  fluid  bottom  in  the  upslope  direction  using  the  incident  plane  wave  reflection 
coefficient  stated  by  Rayleigh.  The  pressure  distribution  at  the  water-bottom  interface  can 
be  used  as  a  source  to  project  beams  into  the  bottom  to  calculate  the  propagation  through 
the  medium.  Brekhovskikh  and  Lysenov  [Ref.  20]  stated  that  for  very  distant  source  and 
image  combinations  or  the  presence  of  many  reflections  allow  image  theory  to  be  applied 
with  Rayleigh  reflection  coefficients  with  insignificant  error.  Baek  [Ref.  21]  then 
developed  a  computer  model  of  pressure  throughout  the  water  column  over  a  fast  botto 
This  model  was  limited  to  a  fast  bottom  due  to  the  critical  grazing  angle  of  the  sound  to 
bottom.  LeSesne  [Ref.  22]  used  the  same  model  as  a  basis  for  a  three-dimensional  model 
which  produced  a  pressure  field  in  the  cross  slope  case  which  was  validated  against 
experimental  results  from  a  model  tank  by  Kosnick  [Ref.  23].  Kaswandi  [Ref.  24] 
developed  a  model  for  a  pressure  field  with  a  slow  bottom  in  downslope.  Nassopoulos 
[Ref.  25]  took  all  previous  models  developed  at  the  Naval  Postgraduate  School  and 
combined  them  into  a  computer  program  called  URTEXT.  He  also  implemented  the  use 
of  the  doublet  expansion  of  the  acoustic  field  in  a  wedge  shqwd  ocean.  Livingood  [Ref. 

26]  then  expanded  the  image  doublets  into  the  cross-slqpe  case.  Joyce  [Ref.  27]  then 
developed  the  quadruplet  expansion  of  the  acoustic  pressure  field  in  the  wedge  shaped 
ocean.  The  quadruplet  expansion  is  a  much  faster  computationally  but  the  approximations 
in  the  reflection  coefficient  and  other  terms  in  the  quadruplet  expansion  equation  produces 
large  errors  at  large  wedge  (P>3')  and  source  angles  due  to  the  angular  differences  between 
the  images  in  the  upper  and  lower  doublets. 
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n.  IMAGE  THEORY 


This  section  is  a  brief  review  of  the  ideas  and  basic  equations  of  image  theory.  The 
following  assumptions  are  made: 

1 .  Sound  velocities  and  densities  are  constant  in  both  the  wedge  and  penetrable  bottom. 

2.  Air-water  boundary  is  a  pressure  release  surface  i.e.  reflections  are  180°  out  of  phase. 

3.  Slope  is  constant  and  all  boundaries  are  planar. 

When  sound  is  transmitted  within  a  waveguide,  there  are  a  number  of  interactions  of  the 
acoustic  waves  with  the  boundaries.  In  image  theory,  each  ray  path  is  replaced  by  an  image 
of  the  source.  The  distance  of  the  source  from  the  receiver  is  the  total  length  of  the  acoustic 
ray.  In  the  case  of  a  sloping  bottom,  the  image  is  placed  equidistant  from  the  reflecting 
boundary  perpendicular  from  the  source.  See  Figure  2-1. 


Figure  2*1  -Wedge  Gemnetry 

For  a  given  bottom  slope  angle  the  number  of  images  N  in  each  of  the  upper  and 
lower  half  spaces  is  given  by 
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-Kt) 


Each  image  is  numbered  from  the  source.  The  source  is  image  number  1  for  the  upper 
half  space.  The  first  image  below  the  bottom  is  image  number  1  for  lower  half  space.  This  is 
continued  until  the  opposite  surface  is  reached.  See  Figure  2-2. 


Figure  l-l  •  Image  Structure  for  Wedge  Shaped  Duct 


The  grazing  angle  6^  of  the  nth  image  sdx>ve  or  below  the  bottom  is  calculated  ficom 

0^  =  np—Y 

a,  *  (n  - 1)^  +  y  for  n  even  (2-3) 


6 


where  y  is  the  angle  of  the  source  from  the  surface. 

The  range  r„  from  each  image  to  the  receiver  is  then  calculated  using 

=  V +  '■2  +  ^0  - cos(fl,  -fi  +  S)  (2-4) 

for  the  upper  series  of  images  and 

r„.  =  +  '■2  +  ^0  - 2rirj cos(d,  +p-S)  (2-5) 

for  the  lower  series  of  images.  In  these  equations,  n  is  the  distance  from  the  wedge  apex  to 
the  source,  r2  is  the  distance  from  the  apex  to  the  receiver,  yo  is  the  cross-slope  or  shoreline 
distance  between  the  source  and  receiver,  and  ^-^+6  is  the  angle  from  the  nth  source  to  the 
receiver  (for  the  upper  image).  See  Figure  2-3. 


Figure  2-3  -  Sooroe-Receiver  Gemnetry  (Upper  Image) 

The  sound  will  f^opagale  from  the  source  to  the  receiver  through  combinations  of  three 
different  path  types:  1)  direct,  2)  surface  reflected,  and  3)  bottom  reflected.  Figure  2-4 
illustrates  a  bottrmi-sutface-bottom  path.  The  upper  drawing  illustrates  the  real  path  of  sound. 


The  surface  is  a  pressure  release  boundary  so  the  incident  angle  is  equal  to  the  reflected  angle. 
The  bottom  has  refractive  [voperties  due  its  conqx>sition  so  the  incident  angle  is  not 
necessarily  equal  to  the  transmitted  angle  but  is  dependent  on  the  p  and  c  ratio  of  the  water 
and  bottmn.  Widi  each  interactitm  with  the  bottom,  a  new  angle  is  produced  which  is  related 
to  die  first  incident  angle.  This  is  discussed  later  widi  the  reflection  coefficient.  The  lower 
drawing  illustrates  the  image  theory  equivalent  The  path  from  the  third  image  to  the  receiver 
is  a  straight  line. 
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The  reflection  coefficient  is  a  {mxluct  of  the  interaction  between  the  boundaries  of  the 


fluid  layers.  The  surface  is  a  pressure  release  boundary  so  each  interaction  produces  a 
reflection  coefficient  of  -1.  Each  interaction  with  the  bottom  requires  calculation  of  the 
reflection  coefficient.  This  coefficient  is  a  function  of  the  speed  of  sound  in  e^h  medium,  the 
density  of  each  medium,  and  the  grazing  angle  of  the  ray  equivalent  on  the  bottom.  The  sine 
of  the  grazing  angle  Omt-  is 


sin 


Q  _  r,  sin(0„  -  2mfi)  +  sin[(2m  - 1)/3  +  ^1 


(2-6) 


for  the  upper  images  and 


r,  sin(fl,  -  Imfi)  +  sin[(2m  + 1))3  - 


(2-7) 


for  the  lower  images.  See  (Figure  2-5). 


n  odd 
n  even 


Fl^iire  2-5  -  bicideiit  Angle  d^m  Calculation  (Upper  Half  Space) 
The  reflection  coefficient  for  each  interaction  R^m  tbcn  is  given  by 
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—sin 


\>l2 


(2-8) 


^sin0„  -ha)  -aj 

where  pi,  P2  are  densities  of  the  water  and  bottom  respectively,  j  =  V^,  and  a  and  b  are 
given  by 


a  = 


\^2J 


-  COS^  0. 


(2-9) 


b  =  2 


r.^  \ 


or 


(2-10) 


For  the  above  equations,  ci,  C2  are  the  speeds  of  sound  of  the  water  and  bottom, 
respectively,  and  (xJk2  is  the  bottom  loss  coefficient. 

The  pressure  for  each  half  space  is  given  by  summing  the  contribution  for  the  upper  and 
lower  images  respectively 

p,=i— (2-11) 

»«i  ^11+  "“I 

where  n«  =  1  when  nsl,2  in  (2-11)  and 

P,  =  i— exp{-Ar..X-l)"^’^'n/?„-  (2-12) 


n»l  ^n- 


msO 


The  number  of  interactions  m  with  the  bottom  is  related  to  the  image  number,  so  the  limit  M 


is  the  integer  of 


'n-H'  _ 
.  2  . 


1  for  the  nth  image.  For  convenience  we  have  omitted  the  factor 


which  would  otherwise  sqrpear  in  (2-11)  and  (2-12). 

Two  special  cases  exist  in  the  upper  half  space.  For  the  case  of  the  direct  path  of  the 
source  to  the  receiver  (nsl),  R  is  1.  In  the  case  of  a  single  reflection  off  the  surface  and  no 
reflections  with  die  bottom  (ns2),  R  is  again  1.  The  phase  inversion  of  the  ray  is  computed  by 
the  (-i)^"'(-«)  term  of  (2-1 1)  and  (2-12). 
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Since  the  model  is  frequency  independent,  a  scaling  distance  Xc  is  used.  Xc  is  defined  as 
the  distance  where  the  lowest  normal  nuxle  attains  cutoff  when  the  bottom  is  fast  (c2>c]). 
This  cutoff  distance  is  calculated  by 


k,X, 


It 

2sin9^tan/3 


ft.  =  cos 


-1 


V^2  > 


(2-13) 


so. 


2k  sin0  iaa$ 

I  c 


(2-14) 


For  a  slow  bottom  an  analogous  scaling,  convenient  for  the  computer  program  is 


k,X,= 


2  tan  ft  tan 


ft  =  cos  ‘ 


(2-15) 


so. 


X  = _ S. _ _ 

c  2k  tan0  taaB 

I  s 


(2-16) 


For  the  total  pressure  at  the  receiver,  the  upper  and  lower  contributions  are  combined, 


thus 


P=P„+P, 


(2-17) 
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in.  DOUBLET  FORMATION 


A.  GENERAL  DESCRIPTION 

An  acoustic  doublet  is  a  pair  of  acoustic  images  180°  out  of  phase.  The  first  doublet 
consists  of  the  source  and  its  reflection,  the  2nd  upper  image.  This  is  called  a  neutral 
doublet  (n=0).  The  reflection  coefficient  is  I  because  its  midpoint  does  not  encounter  a 


Figure  3-1  •  Doublet  Image  Structure 

boundaiy.  The  3rd  and  4th  upper  images  are  the  1st  upper  doublet  (n=:l).  The  Sth  and  6th 
ujqper  images  are  the  2nd  upper  doublet  (n=2).  The  1st  and  2nd  lower  images  are  the  1st 
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lower  doublet  (n=l)  and  the  3rd  and  4th  lower  images  are  the  2nd  lower  doublet  (n=2). 
This  procedure  is  continued  until  the  opposite  surface  is  reached.  See  Figure  3-1. 

The  total  number  of  acoustic  doublet  pairs  is  given  by 


N  =  Int 


Int] 


^360°^ 


P 


-2 


(3-1) 


At  worst,  only  three  images  will  be  missed  and  they  will  be  in  the  higher  order  images  of 
the  source  so  their  contribution  is  negligible  due  to  all  the  reflection  coefficients 
encountered  [Ref.  25]. 


B.  DEVELOPMENT 

The  sound  propagation  for  a  spherically  spreading  point  source  is 

p  =  (3-2) 

r 

where  A  is  the  pressure  amplitude,  r  is  the  range  from  the  source  and  k  is  the  wave 
number  co/c  [Ref.  1 1].  When  two  sources  are  combined  to  form  a  doublet,  the  above 
equation  can  be  expanded  to 

P  _  At.  )  -  Az-  ^  (3-3) 

r.  4 

where  the  subscripts  “+”  and  refer  to  the  upper  and  lower  sources  of  the  doublet  and 
r^  =  r-»-Ar  A^=A-l-M  (3-4) 

r  =  r  —  At  A_  —  A  —  AA 

Tl»  delta  values  are  the  incremental  difference  between  each  source  compared  to  the 
theoretical  point  source.  Figure  3-2  iUustrates  the  doublet  close-up.  The  distance  d  is  the 
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r 

I 

I 

i 


separation  of  the  images  of  the  doublet  pair  and  r  is  the  distance  of  the  doublet  midpoint  to 
the  receiver. 


P(r.  0,t) 


Figure  3-2  •  Doublet  Close-Up 
The  ptessuie  equation  can  now  be  rewritten  as 


rA-l-AA^ 

rA-AA^ 

A 

A 

r- At 

r+Ar 

1  /■  J 

1  r  J 

1.  Balanced  Doublet 


(3-5) 


The  1st  upper  image  has  an  reflection  image  of  equal  amplitude  and  opposite 
phase  which  is  the  second  upper  image.  In  the  far  field,  the  image  pair  is  treated  as  a  single 
doublet  source.  See  Figure  3-3. 
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P  =  2  j—siaikr.  sin(r)sin(0)y  (3-10) 

r 

2.  Unbalanced  Doublet 

For  the  unbalanced  doublet,  the  amplitude  relationship  is  expressed  as 


= 


_  -t" 

2 


AA  = 


A^-A. 

2 


(3-11) 


The  unbalanced  doublet  can  be  viewed  as  a  combination  of  a  balanced  doublet  on  which  is 


superimposed  a  like-charge  doublet  [Ref.  28].  See  Figure  3-4 


The  pressure  equation  is 


P=  2y^sin(kr,sin(y)sin(fl))  +  2^^cos(ifcr,  sin(y)sin(0))  *'’(3-12) 


which  can  be  rewritten  as 
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P  =  2j 


AA 

sin(A;r,  sin(y)sm(0))  -  j — cos(jkr,  sin(y)sin(0)) 

A. 


jica-kr) 


(3-13) 


3.  Distance  firom  Doublet  to  Receiver 

The  distance  from  the  doublet  formation  to  the  receiver  can  be  calculated  using  the 
law  of  cosines.  For  a  three  dimensional  wedge,  the  following  equation  is  used 


rl^  =  rf  -t-  Tj  +  >0  -  2/’,rj  cos(2n/3  ±  8)  (3- 14) 

where  rj  and  r2  are  the  ranges  from  the  apex  to  source  and  receiver  respectively,  yo  is  the 
cross-slope  range,  P  is  the  wedge  angle,  5  is  the  source  angle,  and  n  is  the  doublet 
formation  number.  The and signs  are  reference  to  the  upper  or  lower  images  in 
the  doublet  formation.  See  Figure  3-5. 


Figure  3>5  •  Doublet  Three  Dimensional  Geometry 

Equation  3-13  can  be  rewritten  by  adding  and  subtracting  2r/r2.  This  yields  the 
following  expression 
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=  r,  -  2r,rj  +  r]  +  >„  +  2r,rj  -  2r,r2  cos(2nj3  ±  5) 

By  combining  terms  and  factoring,  the  following  is  the  result 

<=('•,-'•2)“  l+7^^r  +  7^ii^(l-cos(2n^±«)) 
(^1-^2)  In  -'•2) 


^  (n-'-2)  (n-^2) 

4.  Pressure  From  Each  Doublet 


(3-15) 


(3-16) 


(3-17) 


A  generalization  of  Equation  3- 1 3  can  be  made  by  substituting  A„  and  for 


and  AA  where 


A  =A+±; 
^  2 

and  6^  for  6  where 


A4.  =  ^* 
^  2 


(3-18) 


e,,  =  2n/3±S 

so  the  pressure  from  each  doublet  is 


(3-19) 


-2j4i  sin(fer, 
'■ntL 


sin(y)sin(0,i))-  j— ^cos(fcr,ysin(0„4)) 


(3-20) 
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IV.  QUADRUPLET  FORMATION 


A.  QUADRUPLET  GENERAL  DESCRIPTION 

A  quadruplet  is  formed  by  combining  the  upper  and  lower  doublets.  The  visualization 


Figure  4*1  -  Qnadruiriet  IDustratkm 


of  the  quadruplet  is  illustrated  in  Figure  4-1.  The  primary  doublet  (the  source  and  its  first 
upper  image)  are  considered  separately.  The  first  quadruplet  consists  of  the  third  and 
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fouith  upper  image  (first  upper  doublet)  and  the  first  and  second  lower  image  (first  lower 
doublet).  The  second  quadruplet  consists  of  the  fifth  and  sixth  upper  images  (second  upper 
doublet)  and  the  third  and  fourth  lower  images  (second  lower  doublet).  The  process  is 
cmitinued  until  the  opposite  surface  is  reached. 

The  number  of  quadruplets  formed  is 


(4-1) 


At  wm^t,  only  3  images  will  be  missed  and  they  are  highest  order  images  so  their 
contribution  is  negligible  due  to  all  the  reflection  coefficients  encountered. 


B.  QUADRUPLET  DEVELOPMENT 
1.  Balanced  Quadruplet 

The  complementary  doublets  are  180‘  out  of  phase.  The  pressure  of  the 
quadruplet  is  given  by  combining  the  nth  upper  and  lower  doublets  from  Equation.  3-10  of 
the  balanced  douUet  to  frxm 


P^=2j^[  sin(ikr, sin(y)sin(0,^))e"^* 
-sin(*r,  sin(  y)sin(0,_ 


(4-2) 


By  substituting  Equatitm  3-19  for  0^  and  using  the  approximation  that  sin(Y)=Y>  the 


balanced  tpiadrupl^  pressure  equatkm  is 


P,  =  27  sin(itr,  ysinilnp + 

r* 

-sin(Ar,  yunilnfi  -  8y)e'^- 


(4-3) 


The  an^litudes  of  the  upper  and  lower  doublets  are  not  equal  because  of  the 
distance  differential  betweoi  them  as  semi  in  Figure  4-2.  This  leads  to  the  concept  of  the 
unbalanced  quadruplet. 
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Figore  4-2  -  mustration  of  UidMdanced  Quadnqilet  Amplitude 


The  At  terms  are  derived  from  the  comparison  of  the  distance  from  the  primary 
doublet  to  the  midpoint  of  the  upper  and  lower  doublets  in  each  quadruplet  (Figure  4-2). 
This  ensures  that  all  phase  angles  are  calculated  from  a  common  reference  point.  Using 
Equatimi  3-17,  by  setting  yopO,  we  can  make  the  following  approximation, 

- 
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(4-5) 


^±  s  |r,  -  rJ  I  +  . -t  (I  -  cos(2n)3  ±  5)) 

1  lr,-r,) 


Now  Arn  can  be  apiKOximated  using  the  following  relation 


Flgiire  4-3  •  lUiistnition  of  Ar  in  a  Qnadrnplet 

The  cosine  tenn  in  the  above  expression  can  be  expanded  using  a  trigonometric 

identity 


cos(2n^  ±  5)  -  cos(2n^)cos(^)  ±  sin(5)sin(2nj3)  (4-7) 

Small  angle  apfnoximatitMis  can  be  used  for  5  where  cos(5)sl  and  sin(5)s5  so 

cos(2nfi  ±  ^)  =  cos(2n)9)  ±  S  sin(2n^)  (4-8) 
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Equation  4-8  can  be  substituted  into  Equation  4-6  which  can  be  substituted  into 


Equation  4-3  to  yield 

P,  =  sin(ifcr,ysin(2n)9  + 

-sin[kr{ysin{2nfi  -  ] 

For  convenience,  the  following  terms  are  defined 
0  =  Jk/i/t(l  -  cos{2nfi)) 

d  =  kr^S  sin(2nfi)  (4- 10) 


Equation  4-9  now  becomes 

[  sia(ir,rsm(2nfi+S))e-''^ 

-sin(/tr,ysin(2n/3-5))e^'*'  ] 

The  amplitudes  of  the  upper  and  lower  images  of  the  doublets  are  not  equal.  By  separating 
the  amplitudes  and  distributing  it  into  the  main  equation  in  a  process  similarly  done  in 
Equaticm  3-12  and  3-13,  the  unbalanced  quadruplet  equation  is  the  result 


[i4),|^sin(itr, 7sin(2n^  +  5))  -  J~‘-cos(kr,ysm(2nJ3  +  5)) (4-12) 
-B,  j^sin(fe-,ysin(2n/3  -  5))  -  J^cos(kr,rsm(2nJ3  - 


where  and  Bn  are  the  anqilitudes  of  the  upper  and  lower  nth  doublets  at  the  midpoint 
and  AAn  and  ABn  are  the  differences  in  anq)litude  of  the  nth  upper  and  lower  doublet 
images. 


23 


V.  REFLECTION  COEFFICIENTS 


A.  REFLECTION  COEFFICIENTS 

The  reflection  coefficient  is  discussed  in  Chapter  2  with  the  image  theory 
development.  The  reflection  coefficient  from  [Ref.  29]  is 


p,  . 
—sin 


—sin 

A 


_ _ 

^2 


It  can  be  rewritten  as 


(5-1) 


bc- 


sip(A) 


bc+ 


sin(g,) 

sin(0«») 


(5-2) 


where  6Ur  is  the  grazing  angle  of  the  nth  image  going  through  the  mth  plane  and 


c  =  ^  (5-3) 

A  Cl 

Sin(6[:)  can  be  expanded  and  the  reflection  coefficient  can  be  rewritten  as 


be- 


<^"cos*(0^) 


R_  sin(g^) 

yi-g  cos"(o 


bc  +  - 


(5-4) 


The  term  under  the  radical  can  be  manipulated  and  written  in  terms  of  the  critical  angle,  6c , 
firom  Equation  2-13 
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(5-5) 


^l-c^cos(eJ=^l-~J^  =c^sia^0j-sm^ej 
Equation  S-S  can  be  combined  into  Equation  S-4  and  simplified. 


By  using  the  substitution 


^■^sin(a^) 

sin(ej 


the  reflection  coefficient  takes  the  form 


(5-6) 


(5-7) 


(5-8) 


B.  SLOW  BOTTOM  REFLECTION  COEFFICIENTS 

A  slow  bottom  is  when  the  speed  of  sound  in  the  bottom  material  is  slower  than  in  the 
water  (ci>C2).  This  is  an  easier  case  to  handle  because  there  is  no  critical  angle  constraint, 
the  reflection  coefficient  is  real  and  diminishes  quickly  to  zero. 

Equation  5-8  can  be  expanded  to  by  using 

sin(eJ  =  l-^  =  ^^  (5-9) 

resulting  in 
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b+ 


1-1 


c^-n  1 


1- 


(5-10) 


c'  jsin'(0^) 
Reversing  the  order  of  the  terms 


By  pulling  out  the  first  term 


(5-11) 


0-^ 

1-c^  1 

’  1  _^2 

1-c^  1 

11  t-  ^  ‘'inf  ft  ) 

sin(0^)\ 

2_^2 

(5-12) 


Dividing  the  tq)  and  bottom  by  lA)  and  making  an  approximation  with  the  square  root 
results  in 


1-. 


1-c^ 


1 


H- 


1  c" 


hsin(e^)V  21-c 


:sin^(^»,) 


(5-13) 


Using  the  relation  of 
l  +  yf 

Equation  S-13  can  be  approximated  by 

2bc  ^ 


(5-14) 


(5-15) 


Substituting 
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a  = 


2bc 

1-c^ 


(5-16) 


The  final  product  is 

(5-17) 

Each  ray  intersects  the  bottom  2n-l  times.  The  product  of  the  reflection  coefficients 
will  be 


(5-18) 


Expanding  the  equation  yields 


*.=  n«le.l= 


(5-19) 


where  rj  =  £  +  S.  Multiplying  exponentials  with  the  same  base  is  equivalent  to  add  their 
exponents.  The  terms  in  can  be  rearranged  to 


By  noticing  that  the  sum  of  the  expontentials  at  a  index  is  equal  to  the  index  squared,  the 
result  is 


(5-21) 

For  small  angles,  the  first  term  in  the  exponential  is  very  small  compared  to  the  second  so 
the  cumulative  reflection  coefficient  for  a  slow  bottom  is 


(5-22) 

C.  FAST  BOTTOM  REFLECTION  COEFFICIENTS 

A  fast  bottom  is  when  the  speed  of  sound  in  the  bottom  is  faster  than  the  speed  of 
sound  of  water  (c2>ci).  At  angles  below  grazing,  the  acoustic  energy  will  propagate  with 
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100%  reflectance  (/?;,;„=1)  until  Oc  (see  Equation  2-13)  is  reached  then  will  decay.  See 
Figure  5-1. 


Figure  5-1  -  Typical  Fast  Bottmn  Cumulative  Reflection  Coefficient  Curve 
The  development  follows  until  Equation  5-8.  By  reversing  the  order  of  terms 


'4-> 

X 

b-j^ 

The  above  equation  can  be  sin^ilified  to 


(5-23) 


1 


-1 


Using  the  Maclaurin  series  [Ref.  30]  of 


(5-24) 
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(5-25) 


In  — =  2tanh-'(z) 

1-z 

Equation  S-22  results  in 

By  taking  out  the  l/x^  term  and  inverting  the  tanh'^ 


[ )t-2taiih"'f y  - ] 


Simplifying  the  above  equation  results  in 


-2taiih" 


or 


-2ytan' 


'te) 


(5-26) 


(5-27) 


(5-28) 


(5-29) 


The  cumulative  reflection  coefficient  is  the  cumulative  product  of  all  previous  reflection 
coefficients  and  will  take  the  form  of 


N  N  -2ycan-'f-;^=] 

R-=n^=n-- 


(5-30) 
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VI.  GRAZING  AND  RECEIVER  ANGLES  IN  QUADRUPLET 


A.  RECEIVER  ANGLE 

The  receiver  angle  is  calculated  using  the  law  of  sines.  The  relation  can  be  seen  in 
Figure  6-1 


li _ _ 


sin(e„i)  sin(2n/3  ±  5) 

By  solving  for  sin(£n)  and  using  the  small  angle  approximation 

e„^  =  -!^sm{2np±S)  (6-2) 

^H± 

Likewise  the  receiver  angle  for  each  image  of  the  doublets  can  be  solved  for.  The  receiver 
angle  for  the  upper  half-space  is  given  by 

_ '[l _ a _ _  (6-3) 

sin(««+±)  sin(2n^  +  5±y) 

sin(6,+4 )  =  -^sin(2nj8  +  5  ±  y)  (6-4) 

^H+t 

Using  a  the  sine  expansion  of  the  sum  of  two  angles,  the  above  equation  results  in 
sin(£,+i)  =  -^[sin(2nj3  +  5)cos(y)±cos(2nj8  +  5)sin(y)]  (6-5) 

^n± 

Making  the  assumption  cos(Y)=l,  sin(Y)=y,  sin(e)=e,  and  substituting  sin(en) 

from  Equation  6-2 
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nth  Bottom 


Figure  6*1  Grazing  and  Receiver  Angles  for  Upper  Doublet 


(6-6) 


en*±  =  €,±-^^ycos(2/i/3-t-<S) 

^rt+ 

for  the  upper  half  space  and 

e«-±  =  e„  ±  —  rcos(2n/3  -  5)  (6-7) 

^11- 

for  the  lower  half  space. 

B.  ANGLE  OF  INCIDENCE  FOR  UPPER  AND  LOWER  DOUBLET 
IMAGES 

The  angle  ftom  the  surface  to  each  virtual  bottom  is 

e^={2n-l)P  (6-8) 

The  angle  of  incidence  of  the  doublet  midpoint  is 

e^^={2n-l)p±S  +  e  (6-9) 

The  angles  of  incidence  with  the  virtual  bottom  is  different  with  each  upper  and  lower 
doublet  image.  The  angles  of  incidence  for  the  ui^r  half  space  are  given  by 

0,,,=(2n-l)/3  +  5  +  e„,,  (6-10) 

By  substituting  Equation  6-5  into  Equation  6-10  and  making  the  standard  small  angle 
assumptions 

=  (2n  - 1))3  +  5 + -^[sin(2/i)3  +  5)(1)  ±  ycos{2nP  -h  5)]  ('6-11) 

^11+ 

Expanding  the  sine  tom  and  making  the  standard  small  angle  assumptions  again 

=  (2n  -\)P  +  S  +  -^[sin(2n/3Xl) + 5cos(2n/?)  ±  ycos{2np  +  5)](6-12) 

^11+ 

Rearranging  the  terms  result  in 
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(6-13) 


=  (2n  -  l)P  +  —sm{2nP)  +  51  +  ^cos(2n)3) 

^n+  .  ^11+ 


±-^  ycos(2n)3  +  5) 

An  assumption  can  be  made  that  if  r,  «  ,  the  sine  term  in  the  above  equation  is  small 


so  the  it  can  be  neglected.  The  result  of  Equation  6-13  is 
- 

6^^^  =  0^-i-S  1  +  ~cos(2nfi}  ±  —  ycos(2nfi  +  5) 
-  ^«+  .  ^"+ 

The  equivalent  equation  for  the  lower  half  space  is 

®<i-±  =  ®itti  ^  ■^cos(2n^)  ±  ycos(2n^  -  5) 

Expanding  the  cosine  term  results  in 


(6-14) 


(6-15) 


+  •^cos{2nfi) 

±-^  y[cos{2n^)cos(5)  -  sin(2«i3)sin(5)] 


(6-16) 


Using  the  standard  small  angle  qrproximatioQs  and  assuming  the  sine  term  is  smaU 


m  1  +  -^cos(2n)3)  ±  ycos(2n^) 

L  ^«+  J  ^«+ 


(6-17) 


for  the  uf^r  half-space  and 


6^_±  =  0^  +  5^1 + -^cos(2»^)  j  ±  ycos(2nj8) 


(6-18) 


for  the  lower  half-space. 


C.  AMPLITUDES 

The  amplitude  of  each  image  in  the  quadruplet  needs  to  be  calculated  due  to  their 
differences.  The  coefficients  of  Equation  6-17  and  6-18  are  used  where 
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(6-19) 


fl,  =s  1  +  ^cos{2nP) 


and 


fr,  =  ■^cos(2n/3) 


(6-20) 


The  reflection  coefficient  can  now  be  written  for  each  image 


-tUtSt 


expanding  the  above  equation  and  using  the  sii»  af>proximation  for  an  angle 
Rb±=Rb«'““**«*''“^* 

Each  of  the  image  amplitudes  can  now  be  expressed 


(6-21) 


(6-22) 


A.,  *  B,,  = 

Using  Equation  3-1 1  and  the  definition  of 


cosh(x)  = 


e'+e'* 


sinh(x)  = 


e*  -  e“* 


(6-23) 


(6-24) 


where  x=noi>n^  gives 


AA.  =  R,  sinh(naB,5)e"^*  AB,  =  R,  sinh(/ia6;,B)e"“"^ 

A,  =  R,  cosh(nab,B)e"®“’*  ^,  =  Rb  cosh(naB,B)e"“'' 
The  equatitms  combine  to  fcNin 


=  tanh(iiab,7) 


B. 


=  tanh(naB,/) 


(6-25) 


(6-26) 


Now  the  full  quadruplet  equati<»  can  be  writtra  as 
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P  »  +  X— R«  cosh{nab,Yy^‘^'''^e-^ 

■  •I 

e'"‘*“'*'^'^(sm(Jkr,ysin(2/i^  +  5))- jtanh(nafc,y)cos(ikr,ysin(2n^  +  S))) 
_^-»aa,s*jiid  ysin{2nfi  -S))-j  tanh(na&„  y)cos{A:r,  y  sin(2n/J  -  5))) 

(6-27) 

where  =  — sin(/tr,  y5) . 

r 
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Vn.  PROCEDURE  AND  RESULTS 


This  thesis  had  three  objectives: 

1 .  Convert  the  WEDGE/URTEXT  to  MATLAB™  and  test  it  in  all  cases  (slow  and 
fast  bottom,  cross-slope,  and  with  bottomloss). 

2.  Linearize  the  quadniplet  expansion  code  by  Joyce  [Ref.  27] 

3.  Extend  the  quadruplet  expansion  to  the  fast  bottom  case  with  a  approximation  of 
die  Rayleigh  reflectitHi  coefficient 

A.  CONVERT  WEDGE/URTEXT  TO  MATLAB 

The  image  model  baseline  program  is  called  WEDGE.  It  is  written  in  GWB  ASIC  by 
Professtv  A.  B.  Coppens  and  is  tun  on  a  IBM  XT(8088)  deslctop  computer.  WEDGE 
was  translated  into  FORTRAN  by  Nassqxiulos  [Ref.  25].  It  was  renamed  URTEXT  and 
is  on  the  mainframe  omqNiter  at  the  Naval  Postgraduate  School.  Translating  it  to  a 
language  that  is  fast  and  can  be  run  on  a  high  performance  desktc^  conqiuter  or 
workstation  was  necessary  to  make  WEDGE  data  readily  available  for  comparison  to  the 
quadtiqil^  expansion.  MATLAB'^  [Ref.  31]  was  chosen  because  it  is  a  high  level 
scripting  language  which  tuns  on  a  variety  of  computer  operating  systems  including  VAX, 
MS-DOS,  UNIX,  and  Macintosh  and  tiiat  it  is  the  mathematics  i»ogram  that  the  author  is 
most  familiar  with. 

VaUdation  of  the  MATLAB^m  version  of  WEDGE,  caUed  WEDGEMAT,  was 
.-inducted  by  comparing  the  results  of  WEDGE  and  URTEXT  to  those  of 
WEDGEMAT. 

Tables  7-1,  through  7-4  show  the  results  of  the  comparison  of  URTEXT  with 
WEDC^MATtr.  WEDGEMATtr  is  the  a  version  of  WEDGE  that  ends  the  calculation 
process  whm  the  pressure  amplitude  is  less  than  0.00000001.  The  major  enw  at  the 
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surface  is  due  to  the  differeiKe  in  precision  of  the  different  languages.  FORTRAN  uses 
single  precision  variables,  accurate  to  six  decimal  places.  MATLAB™  is  precise  to  13 
places,  therefore  at  the  small  pressures  near  the  surface  the  round-off  errors  are  extreme.  It 
should  be  noted  that  the  results  of  Nassoloulos  are  normalized  pressures  where  the  actual 
pressure  is  multiplied  by  the  source  scaling  distance  ri. 

Tables  7-5  through  7-8  show  the  comparison  of  WEDGE  to  WEDGEMATtr.  The 
WEDGE  program  in  GWBASIC  uses  variables  precise  to  12  decimal  places  so  it 
compares  more  favorably  with  the  MATLAB  version.  Runs  were  done  with  fast  bottoms, 
bottomloss,  and  cross-slope  to  test  accuracy.  Tte  results  are  very  favorable. 

A  run  was  made  on  an  10°  slope  with  the  data  taken  from  the  surface  to  bottom 
sampled  every  0.05°.  Figure  7-1  shows  pressure  falls  off  linearly  as  the  sound  nears  the 
surface  (bottom  of  graph).  This  implies  the  model  is  valid  at  the  surface. 


F^pre  7*1  •  Pressure  vs.  Receiver  Angle 
WEDGEMATtr  has  imdergone  basic  proofing.  When  compared  with  previous 
versions  in  FORTRAN  and  GWBASIC,  it  is  more  flexible  both  in  platforms  run,  model 
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variables,  and  speed.  Most  importantly,  now  it  can  be  used  to  check  the  quadruplet 
expansion  data  to  verify  results. 


TABLE  7-1  SLOW  BOTTOM  COMPARISON  OF  URTEXT  AND 

WEDGEMATTR 


P=3  7^=1  ri=l  r2=10  yo=0  pi/p2=0.9 

ci=1500m/s  c2=1475m/s  a/k2=0 _ 


Pressure  Amplitude 

1 

Phase  Angle  I 

Receiver 

Image 

Absolute 

Image 

Absolute 

Angle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.03* 

0.000123 

0.000126 

0.000003 

-54.7* 

48.8* 

5.9* 

0.30* 

0.001255 

0.001250 

0.000005 

■49.8* 

-48.7* 

1.1* 

0.60* 

0.002415 

0.002400 

0.000015 

-48.8’ 

-48.5* 

0.3* 

0.90‘ 

0.003352 

0.003356 

0.000004 

-49.4* 

48.1* 

1.3* 

1.20* 

0.004028 

0.004039 

0.000011 

-48.5* 

47.5* 

1.0* 

1.50* 

0.004384 

0.004390 

0.000006 

-47.5* 

46.4* 

i.r 

1.80* 

0.004369 

0.004377 

0.000008 

■AST 

44.5* 

1.2* 

2.10* 

0.004010 

0.004002 

0.000008 

-42.r 

41.1* 

1.0* 

2.40* 

0.003339 

0.003316 

0.000(H3 

-35.2* 

-34.5* 

0.7* 

2.70* 

0.002471 

0.002467 

0.000004 

-20.8* 

-19.7* 

1.1* 

Ava  Piressure  Diff= 

0.000008 

Avg  Phase  Diff  - 

1.5* 

Note:  Normalized  pressure  equation  used  where  P  = 


+  P^ 

u 
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TABLE  7-2  SLOW  BOTTOM  COMPARISON  OF  URTEXT  AND 

WEDGEMATTR 


P=4  Y=1  ri=5  r2=12  yo=0  pi/p2=0.7 

ci=1500m/s  c2=1420m/s  a/k2=0 _ 


Pressure  Amplitude 

- 1 

Phase  Angle  II 

Receiver 

Image 

Absolute 

Image 

Absolute  1 

Ansle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.04‘ 

0.040927 

0.041003 

0.000076 

-149.3' 

-149.1' 

0.2' 

0.40* 

0.394236 

0.394324 

0.000088 

-148.6' 

-148.5' 

o.r 

0.80* 

0.701851 

0.701860 

0.000009 

-146.5' 

-146.3' 

0.2' 

1.20* 

0.873027 

0.873049 

0.000022 

-142.3' 

-142.2' 

O.r 

1.60* 

0.915180 

0.915294 

0.000114 

-135.2' 

-135.1' 

o.r 

2.(W 

0.887048 

0.887038 

0.000010 

-124.6* 

-124.5' 

0.1' 

2.40’ 

0.851987 

0.851901 

0.000086 

-112.8' 

-112.7' 

0.1' 

2.80* 

0.810276 

0.810180 

0.000096 

-104.5' 

-104.4' 

or  1 

3.20* 

0.704744 

0.704662 

0.000112 

-103.6' 

-103.5' 

o.r 

3.60* 

0.503485 

0.503407 

0.000078 

-115.2' 

-115.0' 

0.2' 

Avs  Ptessuie  Oiffs 

0.000069 

Avg  Phase  Diff  = 

o.r 

Note:  Normalized  pressure  equation  used  where  P  =  r^^P^ 


TABLE  7-3  SLOW  BOTTOM  COMPARISON  OF  URTEXT  AND 

WEDGEMATTR 


P=5  y=4  ri=2  12=400  yo=0  pi/p2=0.6 

ci=1500m/s  c2=1410tn/s  a/k2=0 


Pressure  Amplitude 

Phase  Angle  | 

Receiver 

Image 

Absolute 

Image 

Absolute 

Angle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.05* 

0.000087 

0.000082 

0.000005 

56.3* 

65.2* 

8.9'  1 

0.50* 

0.000818 

0.000812 

0.000006 

62.6* 

65.3* 

2.7*  1 

1.00* 

0.001532 

0.001549 

0.000017 

62.2* 

65.4* 

3.2*  1 

1.50* 

0.002151 

0.002144 

0.000007 

63.0* 

653* 

2.5'  1 

2.00* 

0.002537 

0.002542 

0.000005 

62.6* 

65.6* 

3.0' 

2.50* 

0.002704 

0.002704 

0 

63.3* 

65.7* 

14* 

3.00* 

0.002639 

0.002613 

0.000026 

63.0* 

65.8* 

bo 

• 

3J0* 

0.002276 

0.002268 

0.000006 

62.4* 

65.9* 

3.5*  1 

4.00* 

0.001699 

0.001690 

0.000009 

63.0* 

65.8* 

2.8*  1 

430* 

0.000926 

0.000916 

0.000010 

60.5* 

65.1* 

5.1* 

Avs  Pressure  Diffs 

0.000009 

Phase  Diffs 

3.r 

Note:  Normalized  pressure  equation  used  where  P  =  fj  +  Pj^ 
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TABLE  7-4  SLOW  BOTTOM  COMPARISON  OF  URTEXT  AND 

WEDGEMATTR 


P=6  7^2  ri=0.8  r2=5  yo=0  pi/p2=0.8 

ci=5l500in/s  c2=1450m/s  ot/k2=0 _ 


Pressure  Amplitude 

Phase  Anele 

Receiver 

I  age 

Absolute 

Image 

Absolute 

Anale 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.06‘ 

0.000372 

0.000377 

0.000005 

-50.2* 

48.3* 

1.9* 

0.60* 

0.003728 

0.003726 

0.000002 

48.9* 

48.6* 

0.3* 

1.20* 

0.007214 

0.007205 

0.000009 

49.9* 

49.8* 

0.1* 

1.80* 

0.010212 

0.010210 

0.000002 

-52.2* 

-52.0* 

0.2* 

2.40* 

0.012559 

0.012557 

0.000002 

-55.3* 

-55.3* 

0* 

3.00* 

0.014134 

0.014136 

0.000002 

-60.2* 

-60.1* 

0.1* 

3.60* 

0.014946 

0.014947 

0.000001 

-67.1* 

-66.0* 

1.1* 

14.20* 

0.015168 

0.015166 

0.000002 

-77.0* 

-76.8* 

0.2* 

4.80* 

0.015224 

0.015224 

0 

-90.5* 

-90.4* 

0.1* 

5.40* 

0.015831 

0.015832 

0.000001 

-107.7* 

-107.6* 

0.1* 

Avk  Pressure  Diifs 

0.000003 

Avk  Phase  Diff  = 

Note;  Nonnalized  pressure  equation  used  where  P  =  /; 
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TABLE  7-5  FAST  BOTTOM  COMPARISON  OF  WEDGE  AND 
WEDGEMATTR  WITH  BOTTOMLOSS 


P=3  7=1  ri=l  r2=10  yo=0  pi/p2=0.9 


ci=1500m/s  c2=1666m/s  oc/k2=0.0(X)l 


Pressure  Amplitude 

Phase  Angle 

Receiver 

Image 

Absolute 

Image 

Absolute 

Angle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.03* 

0.01256 

0.01255 

0.00001 

-3.1* 

-3.1’ 

0 

0.30’ 

0.12379 

0.12377 

0.00002 

-3.3* 

-3.3* 

0 

0.60‘ 

0.23680 

0.23681 

0.00001 

-3.8* 

-3.9* 

0.1* 

0.90* 

0.32981 

0.32982 

000001 

4.9* 

4.9* 

0 

1.20* 

0.39269 

0.39272 

0.00003 

-6.2* 

■6.2* 

0 

1.50* 

0.42213 

0.42214 

0.00001 

-8.9* 

* 

00 

00 

1 

0.1* 

1.80* 

0.40586 

0.40585 

0.00001 

-10.7* 

o 

I 

0.1* 

2.10* 

0.34833 

0.34832 

0.00001 

-16.6* 

-16.7* 

0.1* 

2.40* 

0.24512 

0.24513 

0.00001 

-15.6* 

-15.7* 

0.1* 

2J0r 

0.14819 

0.14816 

0.00003 

-13.5* 

-13J* 

0 

Avk  Pressure  Diff= 

0.00002 

Avg  Phase  Diff  = 

0.5*  1 
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TABLE  7-6  FAST  BOTTOM  COMPARISON  OF  WEDGE  AND 


WEDGEMATTR  WITH  BOTTOMLOSS 

P=3  7^1  ri=l  r2=100ycF0  pl/p2=0-9 

ci=1500m/s  c2=1875m/s  a/k2=0.0001 _ 


Pressure  Amplitude 

Phase  Angle 

Receiver 

Image 

Absolute 

Image 

Absolute 

Angle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.03* 

0.00210 

0.00199 

0.00011 

35.4" 

40.3* 

4.9* 

0.30* 

0.02069 

0.02062 

0.00007 

34.9" 

35.4* 

0.5* 

0.60" 

0.03962 

0.03966 

0.00004 

34.8" 

34.9* 

0.1* 

0.90" 

0.05529 

0.05535 

0.00006 

34.5" 

34.5* 

0 

1.20" 

0.06632 

0.06640 

0.00008 

34.1" 

34.1* 

0 

1.50* 

0.07174 

0.07183 

0.00009 

33.6" 

33.6* 

0 

1.80" 

0.07095 

0.07105 

0.00010 

33.  r 

33.1* 

0 

2.10" 

0.06385 

0.06394 

0.00009 

32.7' 

32.7* 

0 

2.40" 

0.05088 

0.05093 

0.00005 

33.5* 

33.4* 

0.1* 

2.70" 

0.02783 

0.02782 

0.00001 

38.9* 

38.9* 

0  1 

Avr 

Pressure  Diff= 

0.00006 

Phase  Diff  = 

0.6*  1 
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TABLE  7-7  FAST  BOTTOM  COMPARISON  OF  WEDGE  AND 


WEDGEMATTR  IN  CROSS-SLOPE  WITH  BOTTOMLOSS 

P=3  Y=1  ri=l  1:2=10  yo=4  pl/p2=0.9 

ci=1500m/s  C2=1666m/s  a/k2=0.0001 


1 

Pressure  Amplitude 

Phase  Anele 

H  Receiver 

Image 

Absolute 

Image 

Absolute 

1  Ansle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.03* 

0.01024 

0.01094 

0.00070 

^.3* 

^3.2* 

i.r 

0.30* 

0.10162 

0.10227 

0.00065 

-44.2* 

44.3* 

o.r 

0.«)‘ 

0.19699 

0.19716 

0.00015 

44.2* 

0.1* 

0.90* 

0.27908 

0.27931 

0.00022 

■44.5* 

44.6* 

o.r 

1.20* 

0.34418 

0.34418 

0 

-44.8* 

44.8* 

0 

1.50* 

0.39920 

0.39896 

0.00024 

-44.8* 

44.9* 

0.1* 

1.80’ 

0.40001 

0.39954 

0.00047 

-39.2* 

-39.3* 

0.1* 

2.10* 

0.35816 

0.35799 

0.00017 

-37.5* 

-37.5* 

0 

2.40* 

0.27680 

0.27659 

0.00021 

-34.5* 

-34.6* 

0.1* 

2.10’ 

0.17050 

0.17036 

0.00014 

-332* 

-33.2* 

0 

Av8  Pressure  Dif^ 

0.00029 

Phase  Diff  s 

0.2* 
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TABLE  7-8  FAST  BOTTOM  COMPARISON  OF  WEDGE  AND 


WEDGEMATTR  IN  CROSS-SLOPE  WITH  BOTTOMLOSS 

P=s3  y=0-5  ri=l  r2=10  yo=8  pi/p2=0.9 
ci=1500m/s  C2=1666m/s  a/k2=0.0001 


Pressure  Amplitude 

Phase  Angle 

Receiver 

Image 

Absolute 

Image 

Absolute 

Angle 

Theory 

MATLAB 

Difference 

Theory 

MATLAB 

Difference 

0.03‘ 

0.00493 

0.00501 

0.00008 

-33.4* 

-27.4* 

6.0* 

0.30V 

0.04887 

0.04911 

0.00024 

-33.4* 

-33.0* 

0.4* 

0.60* 

0.09466 

0.09495 

0.00029 

-33.5* 

-33.4* 

0.1* 

0.90* 

0.013453 

0.13463 

0.00010 

-33.5* 

-33.6* 

o.r 

1.20* 

0.16917 

0.16949 

0.00032 

-34.7* 

-34.6* 

o.r 

1.50* 

0.18791 

0.18845 

0.00054 

-29.0* 

-28.9* 

0.1* 

1.80* 

0.18490 

0.18554 

0.00064 

-27.8* 

-27.8* 

0 

2.10* 

0.16321 

0.16382 

0.00061 

-24.6* 

-24.6* 

0 

2.40* 

0.12657 

0.12706 

0.00049 

-23.4* 

-23.5* 

0.1* 

2.70* 

0.07764 

0.07792 

0.00028 

-21.8* 

-21.9* 

0.1* 

Avg  Pressure  Diff= 

0.00036 

Avg 

Phase  Diff  = 

o.r 

B.  LINEARIZATION  OF  QUADRUPLET  EXPANSION  PROGRAM 

The  original  quadruplet  expansion  program  was  written  by  Joyce  [Ref.  27].  His 
program  is  used  as  a  basis  for  the  fast  bottom  case  examined  later  in  this  chapter.  The 
linearized  {nrogram  is  in  Appendix  A-1  and  A-2  and  has  incorporated  into  it  the  fast 
bottom,  reflection  coefficients  of  Equations  7-1  and  5-28. 
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The  original  quadruplet  expansion  code  used  "For  Loops"  to  compute  its  products. 
The  code  was  linearized  to  speed  processing  time.  A  speed  factor  of  8  was  thieved  with 
the  same  ouq)ut  results  [Ref.  27].  Table  7-9  shows  relative  execution  time  for  a  3°  wedge. 
This  table  is  only  meant  to  give  the  reader  a  'feel'  of  the  time  involved  in  running  the 
programs.  The  WEDGE  program  in  GWBASIC  was  run  on  a  IBM  XT  (8088)  with 
FPU.  All  other  programs  were  run  on  a  Mac  Ilsi  (68030  @  20  MHz)  with  Applied 
Engineering  32k  RAM  cache  card  (with  20  MHz  FPU)  using  MATLAB™  3.5. 


TABLE  7-9  EX) 

ECUnON  TIME  OF  IMA( 

SE  PROGRAMS 

Bottom  Type 

WEDGE 

(truncated 

version) 

WEDGEMAT 

WEDGEMATtr 

QUAD 

Linearized 

QUAD 

Slow 

16s 

42s 

2s 

8s 

Is 

Fast 

124s 

43s 

5s 

8s 

Is 

Another  method  to  speed  calculation  is  to  use  the  "For  Loop"  and  put  in  a  "Break" 
command  to  stop  the  calculation  process  after  the  pressure  for  a  certain  image  or  image 
doublets,  in  the  quadruplet  case,  falls  below  a  specified  level.  This  was  done  in  the 
WEDGE  (truncated  version)  and  WEDGEMATtr.  This  especially  useful  in  a  slow 
bottom  case  where  the  first  3  of  30  quadruplets  (for  a  3’  wedge)  would  be  98  percent  of  the 
amplitude.  This  idea  was  not  explored  for  the  quadruplet  expansion  case. 

C.  FAST  BOTTOM  REFLECTION  COEFnCIENT  APPROXIMATION 
The  Rayleigh  reflectitm  coefficient  without  absorption  is 
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P2^2 


cos"(^^)-l 


+j- 


"1  y 


sin(e^) 


P2^2 


-;• 


\^\ ) 


cos'(®».)-l 


A^i  ■  sin(d„) 
can  be  approximated  by  Equation  S-28  which  is 


(7-1) 


For  a  fast  bottom  3’  wedge  with  the  following  parameters 

c,  =  15(X)m/s  p^—\Kgln^ 

C2*1666m/j  Pj  =l.lllli:g/m^ 

the  reflection  coefficients  are  displayed  in  Table  7-10.  In  a  3*  wedge,  30  reflection 
coefficients  would  be  calculated.  A  fast  bottom  requires  the  first  10  terms  since  the 
cumulative  product  is  multiplied  with  the  amplitude  of  the  quadruplet.  The  cumulative 
product  after  the  tenth  term  is  veiy  small  and  contributes  little  to  the  final  pressure  as  can  be 
seen  in  Figure  5-1. 
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TABLE  7-10  COMPARISON  OF  RAYLEIGH  AND  ARCTAN 
APPROXIMATION  REFLECTION  COEFFICIENTS 


m 

Rayleigh  Reflection 

Coefficient 

Tan-1 

Approximation 

Amp 

Phase 

Amp 

Phase 

1 

1 

149.7° 

1 

149.3° 

2 

1 

121.2° 

1 

117.8° 

3 

1 

95.1° 

1 

83.7° 

4 

1 

70.9* 

1 

38.3° 

5 

1 

46.1° 

0.380 

0 

6 

0.9084 

0 

0.2471 

0 

7 

0.4839 

0 

0.1885 

0 

8 

0.3750 

0 

0.1568 

0 

9 

0.3120 

0 

0.1375 

0 

10 

0.2669 

0 

0.1251 

0 

The  correlation  of  the  two  reflection  coefficients  are  good  until  the  third  quadruplet. 
Then  the  Tan'^  2q>proximati<m  fall  off  much  mcHe  quickly  than  the  Rayleigh  values.  This 
can  be  seen  in  die  conqilex  plane  plot  of  the  reflection  coefficients  in  Figure  7-2.  The  first 
five  reflecticm  coefficients  are  numbered  in  each  case. 
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Figure  7>2  •  Comparison  ct  Complex  Rrilection  Corifkioits 
TABLE  7-11  COMPARISON  OF  THE  CUMULATIVE  RAYLEIGH  AND 


ACTTA 

kN  APPROXIMATION  REFLE( 

CTION  COEFnCIENTS 

Rayleigh  Reflection  CbefScient 

Tan-1 

Approximation 

Quad  Number 

Amp 

Phase 

Amp 

Phase 

1 

1 

149.7* 

1 

149.3* 

2 

1 

-89.1* 

1 

-92.9’ 

3 

1 

6.1* 

1 

-9.2’ 

4 

1 

77.0* 

1 

29.0* 

5 

1 

123.1* 

0.3880 

29.0’ 

6 

0.9084 

123.1* 

0.0959 

29.0’ 

7 

0.4396 

123.1* 

0.0181 

29.0’ 

8 

0.1648 

123.1* 

0.0028 

29.0’ 
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The  most  inqwitant  factor  is  the  cumulative  rci  tdction  coefficient  which  is  shown  in 
Table  7-11.  There  is  a  wide  disparity  between  the  actual  and  ^roximated  reflection 
coefficients.  Hie  correlation  is  good  until  the  second  quadruplet,  then  the  products  diverge. 
The  reastm  for  this  is  the  inaccuracy  of  each  grazing  with  the  bottom  is  amplified  by  the 
complex  multiplicatitMi  of  the  cumulative  product. 

A  test  was  run  comparing  the  amplitude  and  phase  at  the  receiver.  The  output  is 
displayed  in  Table  7-12.  Data  were  taken  from  the  surface  to  the  bottom  in  0.5° 
increments.  The  differences  are  attributed  to  the  inaccuracy  of  the  reflection  coefficient 
after  the  second  quadruplet.  With  the  fast  bottom,  the  equations  are  very  sensitive  to  the 
phase  of  the  reflection  coefficients.  Even  though  the  reflection  coefficients  have  an 
anqilitude  of  1,  the  phase  can  cause  constructive  and  destructive  interference.  Also  with  a 
fast  bottom,  many  of  the  terms  in  the  quadruplet  expansion  are  conqilex  there  by 
magnifying  the  effect  of  the  phase  inaccuracy. 

TABLE  7-12  COMPARISON  OF  AMPLITCDE  AND  PHASE  OF  USING 

RAYLEIGH  AND  ARCTAN  APPROXIMATION _ 

RsQtleigh  Reflection  Coefficient 

Amp  I%ase 


0* 

0 

0 

0 

0 

0.5* 

0.0082 

173.0* 

0.0025 

-83.3* 

1.0* 

0.0164 

172.9* 

0.0049 

-83.2* 

1.5* 

0.0245 

172.7 

0.0074 

-83.1* 

2.0* 

0.0239 

172.4* 

0.0070 

-82.8 

2.5* 

0.0161 

172.1* 

0.0043 

-82.1° 

Tan'l 

Approximation 
Amp  I%ase 


Receiver  Angle 
(from  Surface) 
5 


3.0' 


0.0010 


162.2 


0.0006 


-65.6' 


Now  one  still  has  to  connate  the  results  of  WEDGEMATtr,  to  the  quadruplet 
expansion.  The  quadruplet  expansion  with  the  Rayleigh  reflection  coefficient  was  tested  in 
the  same  conditions  as  in  the  above  table.  The  results  are  displayed  in  Table  7-13. 


TABLE  7-13  COMPARISON  OF  AMPLITUDE  AND  PHASE  OF 
QUADRUPLET  EXPANSION  USING  RAYLEIGH  AND  IMAGE  MODEL 


Receiver  Angle 

(from  Surface) 

5 

Rayleigh  Reflection  Coefficient 

WEDGEMATtr 

Amp 

Phase 

Amp 

Phase 

0* 

0 

0 

0 

0 

0.5* 

0.0082 

173.0* 

0.0328 

-12.3° 

1.0* 

0.0164 

172.9* 

0.0592 

-12.6° 

1.5* 

0.0245 

172.7 

0.0747 

-12.6° 

o 

o 

0.0239 

172.4* 

0.0720 

-5.7 

2.5* 

0.0161 

172.1* 

0.0425 

-1.5° 

3.0* 

0.0010 

162.2* 

0.0019 

19.0° 

The  comparison  of  the  image  model  (WEDGEMATtr)  and  the  quadruplet  expansion 
using  the  Rayleigh  reflection  coefficient  reveals  an  error  factor  of  about  3  in  amplitude. 
This  can  be  attributed  to  the  many  iqiproximations  made  in  deriving  the  equation  for  the 
quadruplet  expansion  which  was  coveted  in  Ch^ter  in  to  VI  of  this  thesis.  Many  of  the 
approximations  ate  first  (rder  and  ate  not  good  enough  to  get  accurate  results. 

The  quadruplet  expansion  using  the  exponential  ^iproximation  for  a  slow  bottom 
(Ref.  27]  is  limited  to  small  wedge  angle,  the  source  close  to  the  surface,  and  the  receiver 
must  be  in  the  'far  field'.  The  reason  for  the  limitation  to  a  small  wedge  angle  (P~3)  is  to 
minimize  the  errors  in  the  small  angle  approximation  (  cos(y)=1  and  sin(5)=5).  The 
source  needs  to  be  close  to  the  surface  (upper  half  of  the  wedge)  so  the  distance  between 
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the  upper  and  lower  doublet  pairs  are  reasonably  close.  The  'far  fip’d'  requires  r,  »  r, , 
which  is  about  SO  dump  distances  (X^).  The  reason  for  this  is  to  insure  the  receiver  angle 
(e)  is  small  and  the  difference  between  the  upper  and  lower  images  in  the  doublet  are  very 
small.  These  requirements  are  needed  to  insure  that  in  the  slow  bottom  case,  where  the 
reflection  coefficient  is  real,  and  the  cumulative  product  of  the  reflection  coefficient  is  ~0  by 
the  third  quadruplet  (of  30)  in  a  3°  wedge.  In  the  fast  bottom  case,  when  the  terms  are 
complex  and  phase  dependent,  the  cumulative  reflection  coefficient  will  go  to  -0  by  the 
ninth  or  tenth  quadmplet.  in  this  latter  case  there  are  many  more  complex  operations  being 
executed  thereby  magnifying  the  need  for  accurate  approximations  of  the  reflection 
coefficients  and  the  terms  in  the  quadruplet  expression  (Equation  6-27). 


Vin.  CONCLUSION  AND  RECOMMENDATIONS 


Image  theory  is  an  ideal  method  for  calculating  the  transmission  loss  in  a  shallow 
water  (wedge  shaped  ocean)  environment.  It  can  be  used  in  cross-slope,  at  all  frequencies 
and  in  transitional  cut  off  regions  that  are  out  of  bounds  for  normal  mode  theories. 

This  thesis  has  successfully  converted  the  GWBASIC  and  FORTRAN  code  of 
WEDGE  and  URTEXT  to  MATLAB™,  which  is  a  high  level  scripting  language.  The 
original  WEDGE  program  is  accurate  in  fast  and  slow  bottoms,  cross-slope,  and  with 
bottomloss.  The  WEDGEMAT  program  retains  all  the  same  qualities  but  is  much  faster 
to  run.  The  quadruplet  expansion  program  for  a  slow  bottom  by  Joyce  [Ref.  27]  was 
linearized  to  increase  processing  speed  by  a  factor  of  8  while  retaining  the  same  results. 

This  thesis  also  incoiporated  a  Tan*^  approximation  of  the  reflection  coefficient  for  a 
fast  bottom  into  the  quadruplet  expansion.  Due  to  the  inaccuracy  of  the  reflection 
coefficient  after  the  second  quadruplet,  the  results  were  not  favorable.  It  was  also 
discovered  that  even  with  an  accurate  Rayleigh  reflection  coefficient,  the  first  order 
approximations  made  in  developing  the  quadruplet  expansion  equation  (Equation  6-27)  are 
not  accurate  enough  for  the  fast  bottom  case. 

The  next  steps  in  this  area  of  research  is  to  find  a  better  approximation  for  the  fast 
bottom  reflection  coefficient  and  to  increase  the  accuracy  of  the  terms  of  the  quadruplet 
expansion  equation  (Equation  6-27). 
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APPENDIX  A-1 


% 

%  QUADATAN 

%  Linearized  Quadruplet  Expansion  of  Method 
%  of  Images.  Based  on  PQUAD.m  by  Joyce. 

% 

%  This  program  will  provide  pressure  amplitude  and 
%  the  phase  angle  at  the  reciever.  All  distances  are 
%  ratios  of  dump  distauices  Xc.  All  cingles  are  inputed 
%  in  degrees.  This  model  includes  Rayleigh  Reflection 
%  coefficient  for  fast  bottom  case  and  exponential 
%  approximation  for  slow  bottom  case. 

% 

%  By  Pat  Takamiya 
% 

clear;  clc;  clg; 

%  Input  Parameters 

B= input (' Enter  bottom  wedge  angle  in  deg  (beta)  ='); 

G= input (' Enter  source  angle  from  surface  in  deg  (gamma)  ='); 
D= input (' Enter  receiver  angle  from  surface  in  deg  (delta) ='); 
dl = input (' Enter  Density  ratio,  water  to  bottom  =’); 

CC=input (' Enter  sound  speed  ratio,  water  to  bottom  =’); 
rl=input (' Enter  range  of  source  from  apex  (rl)  ='); 
r2=input( 'Enter  range  of  receiver  from  apex  (r2)  ='); 
y0=0; 

% 

%  Determine  Number  of  Image  Quadruplets 
%  Convert  Angles  from  Deg  to  Rads 
Nl=fix{90/B)  ; 

B=B* (pi/180) ; 

G=G*(pi/180) ; 

D=D* (pi/180) ; 

% 

%  Determine  Scaling  for  Fast  or  Slow  Bottoms 
tgg=tan(B) ; 

if  CC<1  %  Fast  Bottom 

tgc[l=acos  (CC)  ; 
tgq2=sin(tqql) ; 
else 

tqql=acos (1/CC) ;  %  Slow  Bottom 

tqq2 = tan ( tqql ) ; 

end 

tqq3 =2  *  tqq2  *  tqq ; 
t4=pi/tqq3; 


%  Scaled  wave  number 


% 

%  Calculate  Constants  and  Primary  Doublet 
% 

alpha=2*{l/dl) /{sqrt{ (CC^2)-1) )  ; 

r=sqrt(  (rl''2)  +  (r2"2)  +  (y0''2)-{2*rl*r2*cos{D)  )  )  ; 

r3=abs (r2-rl) ; 

r4=rl*r2/  (r3''2)  ; 

mu=r2/r3 ; 

q=t4*rl*G*D; 

pl= (2/r) *sin(q) ; 

f=0; 

% 

%  Quadruplet  Simnmation 
% 

n=(l:Nl) ; 

s  (n)  =n; 
th(n) =2*n*B; 

d(n)=t4*rl*D*sin(th{n) ) ; 

phi (n) =t4*rl*mu* (1-cos (th{n) ) ) ; 

%  Reflection  Coeff 

ru(n)=sqrt(  (rl^2)  +  (r2^2)  +  (y0'^2)-2*rl*r2*cos(th(n)+D) )  ; 
rl(n)=sqrt{ {rl^2)+(r2^2)+{y0"2)-2*rl*r2*cos(th(n)-D) ) ; 
r5(n)=rl./ru(n) ; 
r6 (n) =rl. /rl (n) ; 
al  (n)  =1  “-rS  (n)  .  *cos  (th(n) ) ; 
bl (n) =r6 (n) . *cos (th(n) ) ; 
if  COl 

rr(n)  =  (B*alpha.*(n.^2) )  ; 

R(n)=-exp(-rr(n) ) ; 
else 

rc=(l/CC)*l/dl; 

ccl=l/CC; 

gl (n) =sqrt (ccl*ccl*cos (th(n) ) -1) ; 
g2  (n)  =sin  ( th  (n) ) ; 
g3  (n)  =  j  *gl  (n)  .  /g2  (n)  ; 

Rl  (n)  =  (rc+g3  (n) )  .  /  (rc-g3  (n) )  ; 

R(n)  =cun®)rod ( -Rl  (n) )  ; 
end 

a(n)=n.*alpha*D.*al(n)  ; 
b(n) =n. *alpha*G. *bl (n) ; 

% 

%  Pressure 

ul (n) =cosh (b (n) ) ; 

u2 (n)=exp(-j*phi (n) ) ; 

u3 (n)=(-2/r) .*R(n) .*ul(n) .*u2(n) ; 

% 
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%  upper  Image 

vl (n) =sin(th(n) +D) ; 
v2(n)=t4*rl*G*vl(n) ; 
v3 (n) =sin(v2 (n) ) ; 
v4  (n)  =cos  (v2  (n) )  ; 
v5  (n)  =-a  (n)  -  j  *mu*d  (n)  ; 
v6  (n)  =exp  (v5  (n) )  ; 
v7  (n)  =tanh  (b  (n) )  ; 

v8  (n)  =v6  (n)  .  *  {v3  (n)  -  j  *v7  (n)  .  *v4  (n)  )  ; 

% 

%  Lower  Image 

wl (n) =sin(th(n) -D) ; 
w2 (n) =t4*rl*G*wl (n) ; 
w3 (n) =sin(w2 (n) ) ; 
w4 (n) =cos (w2  (n) )  ; 
w5  (n)  =a(n)  +j*mu*d(n)  ; 
w6  (n)  =exp  (w5  (n) )  ; 

w7  (n)  =w6  (n)  .  *  (w3  (n)  -j  .  *v7  (n)  .  *w4  (n) )  ; 

% 

%  Summation 

p  (n)  =u3  (n)  .  *  (v8  (n)  -w7  (n) )  ; 
p2=sum(p(n) ) ; 

% 

% 

ps=pl+p2 ; 

phl=atan2 { imag (ps) , real  (ps ) )  ; 
phase= 1 8  0  *phl /pi 
P=sqrt(real{ps)  ^^2  +  imag(ps)''2) 
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APPENDIX  A-2 


% 

%  QUADATAN 

%  Linearized  Quadruplet  Expansion  of  Method 
%  of  Images.  Based  on  PQUAD.m  by  Joyce. 

% 

%  This  program  will  provide  pressure  amplitude  and 
%  the  phase  angle  at  the  reciever.  All  distances  are 
%  ratios  of  dump  distances  Xc.  All  angles  are  inputed 
%  in  degrees.  This  model  includes  ArcTan  Approximation 
%  Reflection  coefficient  for  fast  bottom  case  and 
e:^onential 

%  approximation  for  slow  bottom  case. 

% 

%  By  Pat  Takamiya 
% 

clear;  clc;  clg; 

%  Input  Parameters 

B= input ( 'Enter  bottom  wedge  angle  in  deg  (beta)  ='); 

G=input (' Enter  source  angle  from  surface  in  deg  (gamma)  ='); 
D= input (' Enter  receiver  angle  from  surface  in  deg  (delta) 

= ' ) ; 

dl= input (' Enter  Density  ratio,  water  to  bottom  ='); 

CC=input (’ Enter  sound  speed  ratio,  water  to  bottom  ='); 
r 1= input ( 'Enter  range  of  source  from  apex  (rl)  ='); 
r2=input ( 'Enter  range  of  receiver  from  apex  (r2)  ='); 
y0=0; 

% 

%  Determine  Nximber  of  Image  Quadruplets 
%  Convert  Angles  from  Deg  to  Rads 
Nl=fix(90/B) ; 

B=B* (pi/180) ; 

G=G* (pi/180) ; 

D=D* (pi/ 180) ; 

% 

%  Determine  Scaling  for  Fast  or  Slow  Bottoms 
tqq=tan(B) ; 

if  CC<1  %  Fast  Bottom 

tqql=acos (CC); 
tqq2=sin(tqql) ; 
else 

tqql=acos (1/CC) ;  %  Slow  Bottom 

tqq2 = tan ( tqql ) ; 

end 
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tqq3  =2* tqq2  *  tqq ; 

t4=pi/tqq3;  %  Scaled  wave  niimber 

% 

%  Calculate  Constants  and  Primary  Doublet 
% 

alpha=2*  (1/dl)  /  (sqrt  ( (CC''2}  -1) )  ; 

r=sqrt( (rl^2)+(r2^2)+(y0"2)-{2*rl*r2*cos(D) ) ) ; 

r3=abs (r2-rl) ; 

r4=rl*r2/  {r3''2)  ; 

mu=r2/r3 ; 

q=t4*rl*G*D; 

pl= (2/r) *sin(q) ; 

f=0; 

% 

%  Quadruplet  Summation 
% 

n={l;Nl) ; 

s(n)=n; 
th(n) =2*n*B; 

d(n) =t4*rl*D*sin{th{n) ) ; 

phi (n) =t4*rl*mu* (1-cos (th(n) ) ) ; 

%  Reflection  Coeff 

ru(n)=sqrt(  (rl^2)  +  (r2''2)  +  (y0^2)-2*rl*r2*cos(th(n)+D) )  ; 

rl(n)=sqrt  ( (rl''2)  +  (r2''2)  +  (y0''2)  -2*rl*r2*cos{th(n)  -D) )  ; 

r5 (n)  =rl. /ru(n) ; 

r6 (n) =rl . /rl (n) ; 

al(n) =l+r5 (n) . *cos ( th (n) ) ; 

bl (n) =r6 (n) . *cos (th(n) ) ; 

if  col 

rr  (n)  =  (B*alpha.  *  (n.''2)); 

R(n)=-exp(-rr(n) )  ; 
else 

x(n) =sin(th(n) ) /sin(acos (CC) ) ; 

Rl{n)=-exp(-2*j*atan(  ( (1/dl)  .*x(n) )  . /sqrt  (l-x(n)  .''2) ) )  ; 
R  (n)  =cunprod  (Rl  (n) )  ; 
end 

a(n) =n. *alpha*D. *al (n) ; 
b(n) =n. *alpha*G. *bl (n) ; 

% 

%  Pressure 

ul (n) =cosh (b (n) ) ; 
u2(n)=exp(-j*phi(n) ) ; 
u3 (n)=(-2/r) .*R(n) .*ul(n) .*u2 (n) ; 

% 

%  Upper  Image 

vl (n) =sin(th(n) +D) ; 
v2 (n) =t4*rl*G*vl (n) ; 
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v3 (n)=sin(v2 (n) ) ; 
v4  (n)  =cos  (v2  (n) )  ; 
v5{n)=-a(n)-  j*mu*d(n)  ; 
v6  (n)  =exp{v5  (n) )  ; 
v7  (n)  =tanh  (b  (n) )  ; 

-8  (n)  =v6  (n)  .  *  (v3  (n)  -  j  *v7  (n)  .  *v4  (n) )  ; 

% 

%  Lower  Image 

wl (n) =sin(th(n) -D) ; 
w2 (n) =t4*rl*G*wl (n) ; 
w3 (n) =sin(w2 (n) ) ; 
w4  (n)  =cos  (w2  (n) )  ; 

w5(n)=a(n)+j*mu*d(n) ; 

w6 (n) =exp (w5 (n) ) ; 

w7  (n)  =w6  (n)  .  *  (w3  (n)  -  j  .  *v7  (n)  .  *w4  (n) )  ; 

% 

%  Summation 

p  (n)  =u3  (n)  .  *  ( v8  (n)  -w7  (n) )  ; 
p2=sum(p(n) ) ; 

% 

% 

ps=pl+p2 ; 

phl=atan2 ( imag (ps ) , real (ps ) ) ; 
phase=180  *phl /pi 
P=sqrt(real(ps)^2  +  imag(ps)^2) 
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APPENDIX  B-1 


% 

% 

%  WEDGE  for  MATIAB 
% 

%  Adapted  from  the  BASIC  version  of  WEDGE 
%  Prof  A.B.  Coppens  and  URTEXT  by  G.  Nassopolous 
%  pz  Not  Normalized 
% 

%  This  program  will  provide  pressure  amplitude  and 
%  the  phase  angle  at  the  reciever.  All  distances  are 
%  ratios  of  dun?>  distances  Xc.  All  angles  are  inputed 
%  in  degrees.  This  model  includes  cross-slope  and 
%  bottomless  terms. 

% 

%  By  Pat  Takamiya 
% 

clear;  clc;  clg; 

% 

B=  input (' Enter  bottom  wedge  cuigle  {beta)=  '); 

G= input (' Enter  source  angle  from  the  surface  (gamma) =  ’); 
D=input( 'Enter  receiver  angle  from  the  surface  (delta) =  ’); 
dl = input ( 'Enter  Density  ratio,  water  to  bottom=  '); 

CC=input (' Enter  speed  ratio,  water  to  bottoms  >); 
rl=input ( 'Enter  range  of  source  from  apex  (rl)=  ' ) ; 
r2sinput ( 'Enter  range  of  receiver  from  apex  (r2)=  ' ) ; 
BL=input (' Enter  bottom  loss  coefficient  (alpha/k2)=  ' ) ; 
yO=input ( 'Enter  cross-slope  range  (Yo)=  '); 

% 

%  Determine  Number  of  Image  Pairs 
%  Convert  Angles  from  Deg  to  Rad 
% 

Nl=fix(180/B) ; 

B=B*pi/180; 

G=B-G*pi/180; 

D=B-D*pi/180; 

% 

%  Determine  Scaling  for  Fast  or  Slow  Bottom 
% 

tqq=tan(B) ; 

if  CC<1,  %  Fast  Bottom 

tqql=acos (CC); 
tqq2=sin(tqql) ; 
else 
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%  Slow  Bottom 


tqql=acos (1/CC) ; 
tqq2 = tan ( tqql ) ; 
end; 

% 

%  Determine  Constants 
% 

tqq3 =2  *  tqq2  *  tqq ; 
t4=pi/tqq3; 
c2=CC^2; 

d2=(rl^2)  +  (r2^2)  +  (y0''2)  ; 
r3=2*rl*r2; 
ql=l/sqrt (2) ; 

% 

%  Define  Lengths  of  arrays 
i=l;fix(Nl/2) ; 
c(i)=zeros(l:fix(Nl/2) ) ; 
e(i)=2eros{l:fix(Nl/2) ) ; 
f (i) =2eros (1: fix(Nl/2) ) ; 

% 

%  Determine  Reuige  to  Receiver  for  each  Image 
% 

n=l:Nl; 

nl=2:2:Nl; 

n2=l;2;Nl; 

tl(nl)=(nl*B)-G; 

tl(n2)=((n2-l)*B)+G; 

r8 (n) =sqrt (d2-r3*cos (tl (n) -D) ) ;  %  3-D  Law  of  Cosines 

r9{n)=sqrt(d2-r3*cos(tl(n)+D) ) ;  %  3-D  Law  of  Cosines 

% 

pl=0;p2=0; 

for  n=l:l:Nl, 

s2=  (-!)''( fix  (n/2) )  ; 


%  Calculate  Reflection  Coefficient  for  Upper  Space 
% 

wl=2*c2*BL; 
il=fix{ (n-l)/2) ; 
for  i=l:l:il, 

s(i)=(rl*sin(tl(n) -2*i*B)+r2*sin(2*i*B-D) ) /r8{n) ; 
if  s(i)>=l, 
s(i)=l; 

end; 

c(i)=sqrt(1.0-s(i)''2)  ; 
t=s(i) /dl; 
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dP  dP  dP  dP  dP  dP 


w0=-c2+c  (i)  ''2  ; 
y=sqrt  ( (w0''2)  +  (wl''2) )  ; 
z=(wO)  ; 

if  y<=z, 
y=z; 
end; 

yl=ql*sqrt (y+wO ) ; 
y2=-ql*sqrt (y-wO) ; 
zl=t-y2; 
z2— "yl / 

z3=zl/ ( (zl''2)  +  (z2'^2)  )  ; 
z4=-z2/{ (zl^2)+(z2^2)) ; 
zl=t+y2; 
z2=yl; 

z5= (zl*z3) - (z2*z4) ; 
z6=(zl*z4)+(z2*z3) ; 
e(i) =z5; 
f  (i)=z6; 
end; 

Determine  Pressure  Contribution  for  the  Upper  Images 


zl=0; z2=0; z3=0; z4=0; z5=l; z6=0; 
if  n>2, 

for  i=l;l;il, 
zl=e(i) ; 
z2=f(i); 
z3=z5; 
z4=:z6; 

z5=(zl*z3)-(z2*z4) ; 
z6=(zl*z4)+{z2*z3); 
end; 

end; 

zl=z5; 

z2=z6; 

t=t4*r8 (n) ; 

z3=cos (t) ; 

z4=-sin(t) ; 

z5=(zl*z3)-(z2*z4) ; 

z6=(zl*z4)+(z2*z3) ; 

pl=pl+ (s2*z5/r8 (n) ) ; 

p2=p2+ (s2*z6/r8 (n) ) ; 

il=il+l; 

Determine  Reflection  Coefficient  for  the  Lower  Images 
for  i=l:l:il, 
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dP  dP  (N> 


s (i) = (rl*sin(tl (n) -2* (il) *B) +r2*sin{2* (il) *B+D) ) /r9 (n) ; 
if  s(i)>l, 
s(i)=l; 
end; 

c(i)=sqrt(1.0-s(i)^2) ; 
t=s (i) /dl; 
w0=-c2+c  (i)  ''2; 
y=sqrt  ( (w0''2)  +  (wl''2) )  ; 
z=  (wO)  ; 

if  y<=2, 
y=z; 

end; 

yl=ql*sqrt (y+wO) ; 
y2=-ql*sqrt (y-wO) ; 
zl=t-y2 ; 
z2=-yl; 

z3=zl/(  (2l"'2)  +  (z2"2)  )  ; 
z4=-22/(  (zl''2)  +  (z2^2)  )  ; 

2l=t+y2 ; 

22=yl; 

25=(2l*z3)-(22*z4) ; 
z6=(2l*z4)+(z2*z3) ; 
e(i)=z5; 
f  (i)=z6; 
end; 

Find  Pressure  Contribution  for  Lower  Images 


zl=0 ; z2=0 ; z3=0 ; z4=0 ; z5=l ; z6=0 ; 
for  i=l:l:il, 
zl=e(i) ; 
z2=f (i) ; 
z3=25; 
z4=z6; 

z5=(zl*z3)-(z2*z4) ; 
z6=(z2*z3)+(zl*z4) ; 
end; 

2l=z5; 
z2=z6; 
t=t4*r9 (n) ; 
z3=cos(t) ; 
z4=-sin(t) ; 
z5= (zl*z3) - (z2*z4) ; 
z6= (z2*z3) + (zl*z4) ; 
pl=pl+ (s2*z5/r9 (n) ) ; 
p2=p2+ (s2*z6/r9 (n) ) ; 

end; 
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»  »  » 


Total  Pressure  at  the  Rec lever 

t5=sqrt  (pl''2+p2''2)  %  Pressure  at  Reciever 

%  Phase  Angle  at  the  Reciever 

phaseang=ateui2  (p2,pl)  *180/pi  %  Phase  Angle 


APPENDKB-2 


% 

%  WEDGE  for  MATLAB 
% 

%  Adapted  from  the  BASIC  version  of  WEDGE 
%  by  Prof  A.B.  Coppens  2Uid  URTEXT  by  G.  Nassopolous 
%  pz  Not  Normalized 
% 

%  This  program  will  provide  pressure  amplitude  and 
%  the  phase  angle  at  the  reciever.  All  distances  are 
%  ratios  of  duunnp  distances  Xc.  All  angles  are  inputed 
%  in  degrees.  This  model  includes  cross-slope, 

%  bottomless,  and  is  triancated  at  pz<0 . 00000001 . 

% 

%  By  Pat  Takamiya 
% 

clear;  clc;  clg; 

% 

B=input (' Enter  bottom  wedge  angle  {beta)=  ’); 

G=input ( 'Enter  source  angle  from  the  surface  (gamma) =  ') 
D=input ( 'Enter  receiver  angle  from  the  surface  (delta) = 
dl= input (' Enter  Density  ratio,  water  to  bottom=  ’ ) ; 

CC= input {' Enter  speed  ratio,  water  to  bottom=  ' ) ; 
rl=input (' Enter  range  of  source  from  apex  (rl)=  ’ ) ; 
r2 = input (' Enter  range  of  receiver  from  apex  (r2)=  ' ) ; 
BL=input (' Enter  bottom  loss  coefficient  (alpha/k2)=  '); 
yO = input ( 'Enter  cross-slope  range  (Yo)=  ' ) ; 

% 

%  Determine  Number  of  Image  Pairs 
%  Convert  Angles  from  Deg  to  Rad 
% 

Nl=fix(180/B) ; 

B=B*pi/180; 

G=B-G*pi/180; 

D=B-D*pi/180; 

% 

%  Determine  Scaling  for  Fast  or  Slow  Bottom 
% 

tqq=tan(B) ; 

if  CC<1,  %  Fast  Bottom 

tqgl=acos (CC) ; 
tqq2=sin(tqgl) ; 
else 

tqql=acos (1/CC) ;  %  Slow  Bottom 
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tqg2=t2Ui(tqql)  ; 
end; 

% 

%  Determine  Constants 
% 

tqq3  =2  *  tqq2  *  tqq ; 
t4=pi/tqq3 ; 
c2=CC^2; 

d2=(rl^2)  +  (r2"‘2)  +  {y0^2)  ; 

r3=2*rl*r2; 

ql=l/sqrt (2) ; 

qa=l;  %  Truncation  Tester 

% 

%  Define  Lengths  of  arrays 
i=l:fix(Nl/2) ; 
c(i)=2eros(l:fix(Nl/2) ) ; 
e(i)=zeros(l:fix(Nl/2) ) ; 
f (i)=2eros(l:fix(Nl/2) ) ; 

% 

%  Determine  Range  to  Receiver  for  each  Image 
% 

n=l:Nl; 
nl=2:2:Nl; 
n2=l : 2 :N1; 

tl(nl)=(nl*B)-G; 

tl{n2)={(n2-l)*B)+G; 

r8 (n) =sqrt {d2-r3*cos (tl (n) -D) ) ;  %  3-D  Law  of  Cosines 

r9{n)=sqrt{d2-r3*cos(tl(n)+D) ) ;  %  3-D  Law  of  Cosines 

% 

pl=0;p2=0; 

for  n=l:l:Nl, 

s2=(-l)''(fix(n/2) )  ; 

% 

%  Calculate  Reflection  Coefficient  for  Upper  Space 
% 

wl=2*c2*BL; 
il=fix((n-l)/2); 
for  i=l:l:il, 

s(i)=(rl*sin(tl(n)-2*i*B)+r2*sin(2*i*B-D) ) /r8(n) ; 
if  s(i)>=l, 
s(i)=l; 

end; 

c(i)=sqrt(1.0-s(i)  ''2)  ; 
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ap  dP  dP  dP  dP 


t=s{i)/dl; 
w0=-c2+c  (i)  ''2; 
y=sqrt  ( {w0‘'2)  +  (wl''2) )  ; 
z=  (wO)  ; 

if  y<=z, 
y=2; 
end; 

yl=ql*sqrt (y+wO ) ; 
y2=-ql*sqrt (y-wO) ; 
zl=t-y2; 
z2=-yl; 

z3=zl/ ( (zl^2)  +  (z2"2)  )  ; 
z4=-z2/( (zl^2)  +  (z2^2) )  ; 
zl=t+y2 ; 
z2=yl; 

z5= (zl*z3) - (z2*z4) ; 
z6= (zl*z4) + (z2*z3) ; 
e(i)=z5; 
f (i) =z6; 
end; 


Determine  Pressure  Contribution  for  the  Upper  Images 


zl=0 ; z2=0 ; z3=0 ; z4=0 ; z5=l ; z6=0 ; 
if  n>2, 

for  i=l:l:il, 
zl=e(i) ; 
z2=f (i); 
z3=z5; 
z4=z6; 

z5=(zl*z3) -(z2*z4) ; 
z6=(zl*z4) +(z2*z3) ; 
end; 

end; 

zl=z5; 

z2=z6; 

t=t4*r8 (n) ; 

z3=cos (t) ; 

z4=-sin(t) ; 

z5=(zl*z3)-(z2*z4) ; 

z6=(zl*z4)+(z2*z3) ; 

pl=spl+  (s2*z5/r8  (n) )  ; 

p2=p2+ (s2*z6/r8 (n) ) ; 

il=il+l; 

Determine  Reflection  Coefficient  for  the  Lower  Images 


67 


dP  dP  dP 


for 

s (i) = (rl*sin(tl (n) 2* (il) *B) +r2*sin(2* (il) *B+D) ) /r9 (n) ; 
if  s(i)>l, 
s{i)=l; 
end; 

c(i)=sqrt(1.0-s(i)^2) ; 
t=s (i) /dl; 
w0=-c2+c  (i)  ''2  ; 
y=sqrt  { (w0''2)  +  (wl'^2) )  ; 
z=  (wO )  ; 

if  y<=z, 
y=2; 

end; 

yl=ql*sqrt (y+wO) ; 
y2=-ql*sqrt (y-wO) ; 
zl=t-y2 ; 
z2=-yl; 

z3=zl/ ( (zl"2)+(z2^2) ) ; 
z4=-z2/(  (zl''2)  +  (z2^2) )  ; 
zl=t+y2 ; 
z2=yl; 

z5= (zl*z3) - (z2*z4) ; 
z6= (zl*z4)  +  (z2*z3) ; 
e(i)=z5; 
f (i)=z6; 
end; 

Find  Pressure  Contribution  for  Lower  Images 


zl=0 ; z2=0 ; z3=0 ; z4=0 ; z5=l ; z6=0 ; 


for  i=l;l:il, 
zl=e(i) ; 
z2=f (i) ; 
z3=z5; 
z4=z6; 

z5= (zl*z3) - (z2*z4) ; 
z6=(z2*z3)+(zl*z4) ; 
end; 


zl=z5; 
z2=:z6; 
t=t4*r9 (n) ; 
z3=cos (t) ; 
z4=-sin(t) ; 
z5= (zl*z3 ) - (z2*z4) ; 
z6= (z2*z3) + (zl*z4) ; 
pl=pl+ (s2*z5/r9 (n) ) ; 
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(»  <»  dP 


p2=p2+(s2*z6/r9(n) ) ; 

%  Truncation  Test 

qa=qa*sqrt( (z5*z5) + (z6*z6) ) ; 
if  qa<0 . 00000001,  break,  end, 

end; 


Total  Pressure  at  the  Reciever 

t5=sqrt  (pl''2+p2''2)  %  Pressure  at  Reciever 

%  Phase  Angle  at  the  Reciever 

phaseang=atan2 (p2,pl) *180/pi  %  Phase  Angle 
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